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ABSTRACT 

The paper describes a new upwind eonservative numerical scheme for special relativis- 
tic resistive magnetohydrodynamics with scalar resistivity. The magnetic field is kept 
approximately divergence free and the divergence of the electric field consistent with 
the electric charge distribution via the method of Generalized Lagrange Multiplier. 
The hyperbolic fluxes are computed using the HLL prescription and the source terms 
are accounted via the time-splitting technique. The results of test simulations show 
that the scheme can handle equally well both resistive current sheets and shock waves 
and thus can be a useful tool for studying phenomena of relativistic astrophysics that 
involve both colliding supersonic flows and magnetic reconnection. 
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1 INTRODUCTION 

In many phenomena of relativistic astrophysics, such as 
AGN, GRBs, quasars, radio galaxies, micro-quasars, pulsars 
and magnetars, compact X-ray binaries etc., the magnetic 
field is a key dynamic component. On one hand, the 
magnetic field drives, accelerates and partially collimates 
relativistic outflows from astrophysical black holes, neutron 
stars, and their accretion disks. On the other hand, magnetic 
reconnection and dissipation is responsible for bright ther- 
mal and non-thermal emission from these flows. Recent years 
have seen a remarkable progress in numerical methods for 
ideal relativistic magnetohydrodynamics fKomissarov 1999' 
Koide et airT999 Komissarov 2001; Koldoba et al. 2002: 
Gainmic ct al. 2003] IDuez et al. "20051 IKoide et al. 1999 
nninos et al. 2006] [Shibata fc Sekuguchi 2005 

2006] 



IDel Zanna et al. 2003 



Anderson et al 

A^it6i7ct'ar2006; Ncilson ct al. 2006' Idizuno ot al. 2007 ; 
Mignonc fc Bode 2006 ^ ^cKinney 2006, Noble et al. 2006. 



Giacomazzo fc Rezzolla 20071 IDel Zanna et al. 2007P and 
many interesting and important simulations have been 
carried out already. Quite often the numerical solutions ex- 
hibited violent magnetic reconnection. Although it is indeed 
very likely to occur in the considered astrophysical phe- 
nomena as the result of non- vanishing physical resistivity of 
plasma, both coUisional and coUisionless, the reconnection 
observed in the simulations is of purely numerical origin. It 
is driven by artificial resistivity arising due to truncation 
errors and hence fully depending on fine details of numerical 
schemes and resolution. A code for resistive RMHD would 
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allow to control magnetic reconnection according to the 
incorporated physical models of resistivity. Moreover, the 
relativistic magnetic reconnection by itself is a sufficiently 
rich and important physical process to warrant the effort 
of developing such a code. The only numerical study of 
relativistic magnetic reconnection so far was carried out by 
Watanabe & Yokovama (| 2006p . However, their paper gives 
no details of their numerical scheme and test simulations 
and therefore it is not clear as to how accurate their results 
are and how robust their numerical method is. 

Since many relevant astrophysical phenomena involve 
shock waves, including the fast magnetic reconnection of 
Petcheck type (ILyubarsky 2005'), a useful code should han- 
dle well not only current sheets and filaments but also shock 
waves. It is well known that codes that do not preserve 
the magnetic field divergence free can become unstable and 
crash in the cases with large spacial gradients. Thus this 
issue must be addressed too. Moreover, in the relativistic 
limit the spacial charge density and the advective current 
can become significant and thus the electric charge conser- 
vation has to be enforced. In this paper we describe the 
results of our efforts to construct a code that satisfies these 
criteria. The equations of resistive RMHD are described in 
Section (2] The equations of the so-called augmented system 
of resistive RMHD, that are designed to handle to enforce 
the divergence free condition for magnetic field and the elec- 
tric charge conservation are presented in Section [S] The rel- 
ativistic Ohm law is explained in Section [l] Section [5] gives 
the details of numerical integration. The test simulations are 
described in Section |6] and our conclusions are summarized 
in Section [T] 
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2 BASIC EQUATIONS 

The covariant Maxwell equations are ( |e.g. Jackson 1979| ) 

*F"^ = 0, (1) 

V/sF'"^ = r, (2) 

where F°'^^ is the Maxwell tensor of the electromagnetic field, 
*F°''^ is the Faraday tensor, and J" is the 4-vector of electric 
current. 

In highly ionized plasma, including pair plasma, the 
electric and magnetic susceptibilities are essentially zero and 
one has 



Fa 



-,al3 



where 



(3) 
(4) 

(5) 



is the Levi-Civita alternating tensor of space-time and eap^jj 
is the four-dimensional Levi-Civita symbol. 

In the coordinate basis of global inertial frame of special 
relativity these equations split into the familiar set 



V-B = 0, 

9tB + V X B = 0, 
V-E = q, 

-dtE + VxB = J 
where 

= F 



2 

j'' = 1 



1 ijk rp 



(6) 
(7) 
(8) 
(9) 

(10) 

(11) 
(12) 



are the electric field, the magnetic field, the electric charge 
density, and the electric current density respectively as mea- 
sured by the inertial observer {eijk = eoijfc is the Levi-Civita 
tensor of space.). These equations are consistent with the 
electric charge conservation 



dtq + V-J ^Q. 



(13) 



In magnetohydrodynamics Maxwell's equations are sup- 
plemented with the equations of motion of matter and the 
continuity equation. In the covariant form the equations of 
motion can be written as 



v.r^" = 

where the total stress-energy momentum tensor, 

^ (m) + (e) 1 

is the sum of the stress-energy momentum tensor of the elec- 
tromagnetic field 



(14) 



(15) 



-(F° 
4^ 



'F 



(16) 



and the stress-energy momentum tensor of matter 
TlZ^^wu^n-'+pg''-'. (17) 



Here where p is the thermodynamic pressure, w(jp,p) is the 
relativistic enthalpy per unit volume as measured in the rest 
frame of fiuid (w includes the rest mass-energy density of 
matter p), and is the fiuid 4- velocity. In the global iner- 
tial frame with time-independent coordinate grid eq |14l splits 
into the energy and momentum conservation laws 



dtP + V-Tl^Q, 

where 



e= -{E'^ + B'^)+w-i'^ 



P 



is the energy density, 

S = ExB A-wy^v, 

is the energy fiux density, 

P = ExB + wy^v 

is the momentum density, and 

n 



(18) 
(19) 

(20) 

(21) 

(22) 



EE - BB + 'W'y'^vv + (-{E'^ + B'^) +p] g (23) 



is the stress tensor. Here 7 is the Lorentz factor, v is the 
velocity as measured by the inertial observer, and g is the 
metric tensor of space. 

The covariant continuity equation is 



Vvpu' = 0, 



(24) 



where p is the rest mass density as measured in the rest 
frame of fiuid. In the inertial frame this reads 



dtp"i + V ■ (p7-y) = 0. 



(25) 



Equations (|6I7I8I9I13I18I19I25P constitute the 3-1-1 PDE 
system of relativistic magnetohydrodynamics in special rel- 
ativity. Once supplemented with equations of state, that re- 
late various thermodynamic parameters of matter, and with 
the Ohm law, that couples matter and the electromagnetic 
field, this system closes. 



3 AUGMENTED SYSTEM 

As well known, the divergence free condition © for the 
magnetic field can be treated as a constraint on the ini- 
tial solution of Cauchy problem because the Faraday equa- 
tion (O will then ensure that the magnetic field remains 
divergence free at any time. The equation of electric charge 
conservation is also not independent and follows from the 
Ampere equation @ and the Gauss law ((S)). These proper- 
ties of the differential equations are not preserved by many 
numerical schemes. Indeed, the most straightforward way 
of constructing a self-consistent finite difference counterpart 
for a differential system like electrodynamics is to ignore all 
constraints (non-evolution equations) and to leave out all 
the supplementary laws like the electric charge conservation 
(otherwise the system of finite difference equations becomes 
over-determined). However, it has been discovered that this 
lack of consistency may lead to strong corruption of numer- 
ical solutions in regions with large truncation errors, like 
strong discontinuities, and even cause code crash. In ideal 
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MHD the divergence free condition has been found particu- 
larly important. A number of techniques has been proposed 
to combat the problem. Here we adopt the so-called Gen- 
eralized Lagrange Multiplier method developed by Munz et 
al. H1999| ). The main idea is to create a new, augmented sys- 
tem of differential equations, that will include only evolution 
equations and will have the same solutions of the Cauchy 
problem as the original system provided the initial solution 
satisfies the differential constraints of the original system. If, 
however, the initial solution does not satisfy the constraints 
then the deviations should decay or at least move away as 
relatively high speed waves. This will ensure that the devia- 
tions caused by truncation errors of a numerical method for 
the augmented system remain small. 

To deal with the divergence free constraint we modify 
eas. (|6l7|l so that they become 

(26) 
(27) 



dt^ + V B = 

dtB + V xE + = 0, 



where $ is a new dynamic variable (pseudo-potential) . From 
these equations it follows that $ satisfies the telegraph equa- 
tion 



C»t $ - Kdt$ + V^-I- = 0. 



(28) 



Thus, $ is transported by hyperbolic waves propagating 
with the speed of light and decays if k > 0. For positive 
K the natural evolution of "I> is toward ${r,t) = (unless 
prevented by boundary conditions) and ea. (|26p shows that 
this final state implies divergence free magnetic field. In fact, 
it is easy to see that the divergence of magnetic field also 
satisfies the same telegraph equation 



c»t (V-B) - Kdt{V-B) + V^(V-B) 



0, 



and thus evolves in the same fashion. 

To deal with the Gauss law we modify eqs. 

read 

- dtE + VxB -V'if ^ J, 



(29) 



9} so they 



(30) 
(31) 



where ^' is another new dynamic variable. From these two 
equations and the electric charge conservation it follows that 
the evolution of ^ is again described by the telegraph equa- 
tion 



(32) 



(Although in principle one could use different constants k for 
$ and ^ this brings no benefit.) Thus, 'I', naturally evolves 
in the same fashion as <I>, ensuring that the electrodynamic 
solution is kept consistent with the Gauss law. Similarly, one 
finds that 

- dt{V-E -q)- KdtiV-E -q)+ \/''{V-E - q) ^ 0. (33) 

Summarizing, the augmented system of relativistic 
MHD is 



a* + V B = 

dtB + VxB + = 0, 

9t* + V-E = g - K*, 

- c»t£; + V X B - V* = J, 



(34) 
(35) 
(36) 
(37) 



dtq + V-J = 0. 

dtpj + V ■ p-yv = 0, 
dte + V-S = 0, 

dtP + v-n = o, 

where 

e = i(S^ + B^) + TO7^ - p 



(38) 
(39) 
(40) 
(41) 



S = ExB + w'yv, 

2 



(42) 
(43) 

P = ExB + wy'^v (44) 
n = EE - BB + w-y'^vv + (^(^^^ + B^) +p) 9 (45) 

Every differential equation of the system is an evolution 
equation and a conservation law (with or without a source 
term) , and there is wealth of numerical methods for such sys- 
tems. For example one could use Godunov's upwind scheme 
(|Godunov 1959|) or one of its numerous higher order "chil- 
dren". This simplicity of numerical implementation is the 
main advantage of the method of Generalized Lagrange Mul- 
tiplier. 

4 OHM'S LAW 

In this paper we consider only the simplest case of scalar re- 
sistivity. In strong magnetic field the resistivity (conductiv- 
ity) becomes anisotropic and the tensor description becomes 
more appropriate. We will consider this case in future. 
The covariant form of scalar Ohm's law is 



(46) 



where a = 1/rj is the conductivity, r; is the resis- 
tivity, and go = —Ivu" is the electric charge density 
as measured in the fiuid frame l|Blackman fc Field 19931 
[Lyutikov fc Uzdensky 2003 [ I. In the general inertial frame 
this reads 

J ^ a^[E + vxB - {E-v)v]+qv, (47) 

whereas in the fluid frame one has the usual Ohm law 

J = cjE. (48) 



In the limit of infinite conductivity (a — > oo) eQ. (l47|l reduces 
to 

E + vxB- {E-v)v = 0. 

Splitting this equation into the components that are normal 
and parallel to the velocity vector one obtains 



E^+vxB = Q 



and 



^ll - {E-v)v = 0. 
These show that E\\ = and thus one has the usual result 
E^ vxB, (49) 



the purely inductive electric field. 

Now consider the reduced form of Ampere's law 

- dtE = J, 



(50) 
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which is of interest for numerical schemes using time- 
splitting technique. When splitted into components normal 
and parallel to the velocity vector this equation reads 



9tB|l + (T7[£;|| - {E-v)v] = 0, 
dtE± + a"/[E± +vxB]=0. 



(51) 
(52) 



The solutions of initial value problem for these linear equa- 
tions are 



E« = Bijcxp ( --t 



and 



E±=El + (B'i - El) cxp i-ajt), 



(53) 



(54) 



where E*j_ = —v x B and sufBx denotes the initial values. 
One can see that for relativistic flows the normal component 
of electric field approaches the inductive value E*^ faster 
than the parallel component approaches zero. 



5 NUMERICAL METHOD 



The evolution equations p4l41|l can be written as conser- 
vation laws. In Cartesian coordinates, and this is the only 
type of coordinates we use in the paper, the system can be 
written as a single phase vector equation 



dQ{V) dT^'iV) 



dt 
where 



dx"" 



(55) 



We have found useful to split the source term into two 



parts 



SaiV) 



—qv 




V 0" / 



and SbiV) 



-K<E> \ 





V 0' J 



where 



a-ylE + vxB - {E-v)v] 



is the conductivity current. The source term Sb is potentially 
stiff (in the case of low resistivity) and is treated via the 
time-step splitting technique by Strang (| 1968p . That is first 
the solution is advanced via integration of equation 



dQjP) 
dt 



■Sb{V), 



(56) 



over the half time-step, Af/2. Then the solution is advanced 
via second-order accurate numerical integration of equation 



dQ(V) ^ dr 



dt 



(57) 



over the full time-step. Finally, the solution is advanced via 
integration of equation (|56p over the half time-step once 
more. Thanks to the fact that all equations in (|56|l are lin- 
ear the integration is carried out using analytic solutions, in 
particular solutions (|53l54p are utilized at this stage. This 
removes stability constraints of the time step otherwise im- 
posed by iSi,. In principle all source terms could be passed to 
eq. (|56|l but this somehow results in reduction of accuracy. 
Equation ([57)l is integrated explicitly 



Q 



q 

PI 

V / 



/ 1 \ 



B' 



S = 



-.r 




0" 



are the vectors of conserved quantities, primitive quantities, 
and sources respectively and 



\ 



e"^''Ek + <S>g'" 
E' 

jm 

pvT 



is the vector of corresponding hyperbolic fluxes. Here = 
7?;' are the spatial components of 4- velocity, g'-* are the com- 
ponents of the metric tensor of space (given by Kronecker's 
delta 5'-' ), e'-'* is the Levi-Civita alternating tensor of space. 



Sn + l = Qn + At 



E 



Ax" 



+ At 5,,„+ 1.(58) 



Here Q„ is the conserved quantity of a cell at t = t„, Qn+i 
is the conserved quantity of this cell nt t = tn + At, i 
is the source term of the cell at t = t„ -f At/2, i „, i is 
the flux though the right interface and J-^_i. is the flux 
though the left interface of the cell normal to the direction 
of x"^ at time t = t„ + At/2. Ax"^ is the ceU size in this 
direction and A^^ is the number of spatial dimensions. 

To determine the sources and fluxes at half time-step 
the solution is temporary advances via 



Q„+ 1 = Sn + 



At 



E 



aJ™ 



+ ^Sa,n. (59) 

The interface fluxes ^^^i „ are computed using the HLL- 
prescription (jHarten et al. 1983^ : 

_ m+TT.n m+TT.n m+7T,n m+TT,ri 
^m+l.n = ^ ^ (60) 
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Figure 1. Stationary fast shock. Left panel: = -fv^; Middle panel: gas pressure (dashed line and stars) and magnetic pressure (solid 
line and circles); Right panel rest mass density, p. 




-0.5 



0.5 



Figure 2. Stationary slow shock. Left panel: = 7U^; Middle panel: gas pressure (dashed line and stars) and magnetic pressure (solid 
line and circles); Right panel rest mass density, p. 



where indexes L and R refer to the states respectively to the 
left and to the right of the interface (which can be considered 
as the location of discontinuity in the solution at t = t„). 
Note the simplification of the general HLL prescription due 
to the fact that the majcimum characteristic speed of the 
system in each direction equals exactly to the speed of light 
(unity in our dimensionless equations). For the auxiliary half 
time-step these left and right states are found via the piece- 
wise constant reconstruction of numerical solution in each 
spatial direction 



Pn^P^ for 



Ax" 



< X < Xc + ■ 



Ax" 



(61) 



where a;™ is the coordinate of the cell center and Pn is the 
phase state vector of the cell. 

The auxiliary solution is then used for another, now 
quadratic reconstruction of numerical solution within each 
cell 



Ax" 



Ax" 



(62) 



Obviously, ai and a2 are the first and the second order 
derivatives of the reconstructed solution and these are to 
be found from the numerical solution using one of many 
existing non-linear limiters (needed to avoid spurious oscil- 
lations). In this particular paper ai is found using the same 
limiter as in our ideal MHD code (|Komissarov 1999|) 



ai = av {P'l,P'r) 
where 



(63) 



_ P^ - Pr + 1 - Pr 

I L - T-Z ' I R - 



av(a, h) 



Ax'" ' ■ " Aa;' 
are the left and right numerical approximations of the first 
derivative (i is the cell index along the direction of a;'"), and 

if a6 < or + 6^ = 0, 
ii^|±f|i if ah>Q and a^' + h'^Q ''^^> 

To find a2 we use a similar procedure. First we compute 
the left, center, and right numerical approximations for the 
second derivative, P'l, Pq, and P'^, and then we feed them 
to the minmod function with three arguments 

a2 = minmod(Pi', Pc,Pr), (65) 
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X 

Figure 5. Self-similar current sheet. The dashed and the dash- 
dotted lines show the exact solution at i = 1 and i = 9 respec- 
tively. The solid line shows the numerical solution at t = 9; it is 
indistinguishable from the exact solution on this plot. 



The electric field is found via the ideal equation 

-E = ~vxB. (68) 

The computational grid is uniform and has 100 cells in 
[—1,-1-1] and the initial solution is set as a discontinuity at 
X = 0. The resistivity is = 0.01. Figure [T] shows the numeri- 
cal solution at t = 3.0 by when the secondary waves created 
during the development of the dissipative shock structure 
have left the grid. One can see that the shock jump is cap- 
tured very well. The fact that there are only 3 grid points 
in the shock structure tells that the shock is unresolved and 
suggests that the shock structure might be dominated by 
numerical dissipation. This is confirmed by the simulations 
with higher resistivity (see fig|3]). 

6.1.2 Stationary Slow Shock 

To set up this test we also solved the ideal relativistic MHD 
shock equations describing stationary shocks. Now the se- 
lected particular solution is 

Left state: 

B = (5.0, 3.511, 0.0), -yv = (0.6082, 0.0, 0.0), p = 1.0, 
p = 10.0, g = 0, $ = 0, * = 0; 



where 



if ab < or 
6c < 



minmod(a,fe,c)=<( ^.^(^,,^^) ^^^^^q (66) 

max{a,b,c) if a,fe, c<0 

The left and right states of each cell interface that are 
found via this second reconstruction are then used to com- 
pute HLL-fluxes ^^^i of equation (|58|) . The resulting 
scheme is second order accuracy in time and third order 
accuracy in space. 



6 TEST SIMULATIONS 

In these test simulations we use the polytropic equation of 
state 



Right state: 

B = (5.0, 2.287, 0.0), -yv = (0.4096, -0.2147, 0.0), p = 1.485, 
p = 17.03; g = 0, $ = 0, * = 0. 

The computational grid is uniform and has 100 cells in 
[—1, -1-1] and the initial solutions is a discontinuity at a; = 0. 
The resistivity is r; = 0.01. Figure [2] shows the numerical 
solution a.t t = 4.0 by when the secondary waves created 
during the development of the dissipative shock structure 
have left the grid. Again, the shock jump is captured very 
well but now there are more then 10 grid point in the shock 
structure. This suggest that shock may be resolved. How- 
ever, the rather small increase in the shock width between 
the cases with r) — 0.005 and rj = 0.01 shows that numer- 
ical dissipation is still important for 77 = 0.01 and only for 
higher resistivity the shock structure becomes fully resolved 
(see figlD). 



p+ jT— j-P (67) 
with the ratio of specific heats F = 4/3. 



6.1 One-dimensional test problems 

6.1.1 Stationary Fast Shock 

To set up this test we solved the ideal relativistic MHD 
shock equations describing stationary shocks. The selected 
particular solution is 

Left state: 

B = (5.0, 15.08, 0.0), jv = (4.925, 0.0, 0.0), p = 1.0, 
p = 10.0, g = 0, $ = 0, * = 0; 

Right state: 

B = (5.0, 28.92, 0.0), -yv = (0.6209, 0.1009, 0.0), p = 7.930, 
p = 274.1; g = 0, $ = 0, * = 0. 



6.1.3 Alfven wave 

To set up this test we utilized the analytical solution for 
ideal MHD Alfven waves obtained in Komissarov (ll997|) . In 
this test p — 1.0, p = 1.0, = 1.0, and the Alfven speed 
Ca = 0.4079. Initially the wave occupies the zone xq < x < 
xi, with Xq — —0.8, xi = 0.0. To the left of the wave B = 
(1.0,0.1,0.0), 71; = 0. In the wave the angle 9 between the 
tangential component of magnetic field and the y-axis varies 
as 



2C'), 



^= {x- xo)/{xi - Xo), 



that gives vanishing first derivatives at xo,i. The initial elec- 
tric field is computed via ea. (l68|l and the electric charge den- 
sity via eq.®. The computational grid is uniform and has 
400 cells in [—1,-1-1]. The resistivity is set to a relatively 
small value, rj — 0.003, in order to get closer to the ideal 
case. The simulations are continued up to t = 1.5 and then 
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Figure 3. Dependence of shock structure on resistivity. Left panel: Stationary fast shock for rj = 0.01 (solid line), 7] = 0.03 (dashed), ?7 = 
0.09 (dash-dotted), r] = 0.18 (dotted), 7] = 0.36 (dash-triple-dotted); Right panel: Stationary fast shock for rj = 0.01 (solid line), r; = 0.02 
(dashcd),r) = 0.04 (dash-dotted), r? = 0.08 (dotted). 




Figure 4. Alfven wave. 



compared with the exact solution of ideal MHD at the same 
time (figSJ. One can see that the agreement is pretty good. 
The ideal solution keeps the wave profile totally invariant, 
however the numerical solution is a little distorted, mainly 
due to numerical dissipation (this is confirmed by studying 
the dependence on rj). 

When the zero gradient boundary conditions (free-fiow) 
are utilised in the simulations then both the fast and the 
slow waves do not get refiected of boundaries and cleanly 
pass through. However, the Alfven waves exhibit noticeable 
refiection (in contract to the results with our ideal MHD 
code). We have not figured out yet as to how to avoid such 
a refiection. 



6.1.4 Self-similar current sheet 

Assume that B — (0.0, B{x, t),0.0), the magnetic pressure is 
much smaller than the gas pressure everywhere, and B{x, 0) 
changes sign within a thin current layer of width Al. Pro- 



vided the initial solution is in equilibrium, p =const, the 
evolution is a slow diffusive expansion of the layer caused 
by the resistivity and described by the archetypal diffusion 
equation 

dtB - r]dlB = 0. 

As the width of the layer becomes much larger than Al the 
expansion becomes self-similar 

B(x,t) =Boerf (^^i=^ , ^^t/x\ (69) 

where erf is the error function, and this analytic results can 
be used to test the resitive part of the code. In the test prob- 
lem that is presented here the initial solution has uniform 
distribution of P = 50.0, p = 1.0, E = 0, and -yv = and 
the initial magnetic field is given by ea. (|69p for Bo = 1-0, 
t — 1, and rj = 0.01. The computational grid is uniform and 
has 200 cells in [—1.5,4-1.5]. The numerical simulations are 
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0,1 o.s 




Figure 6. Strong cylindrical explosion. Top left panel: and magnetic field lines; Top right panel: £?" and magnetic field lines; Bottom 
left panel: logiop, gas pressure; Bottom right panel: Lorentz factor. 




Figure 7. Strong cylindrical explosion. This plot show slices along x = (dashed line and stars) and y = (solid lines and circles) for 
gas pressure (left panel), Lorentz factor (middle panel), and B^ (right panel). 
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continued wp to t — 8 and then the numerical solution is 
compared with the solution (|69|) for t — 9. The results are 
shown in figure [5] - one cannot see the difi'erence between 
the solutions. 

6.2 Multi-dimensional tests 

All the one-dimensional problems, that are described above, 
have been used to test both the 2D and 3D versions of the 
code via application in all two/three directions of the Carte- 
sian grid. The results are almost identical to that of ID 
tests. In addition we considered several generically multi- 
dimensional problems. 

6.2.1 Strong cylindrical explosion 

Strong symmetric explosions are useful standard tests for 
MIfD codes even if there are no exact analytic solutions to 
work with. This is because the generated shocks make all 
possible angles to the grid and to the magnetic field thus 
allowing to detect well hidden bugs and to reveal potential 
weaknesses. In this problem the Cartesian computational 
domain is (-6.0, -1-6.0) x (-6.0, +6.0) with 200 equidistant 
grid points in each direction. The initial explosion zone is a 
cylinder of radius r = 1 centered onto the origin. Its pressure 
and density are set to p = 1 and p = 0.01 for r < 0.8 and 
exponentially decrease for 0.8 < r < 1.0. The ambient gas 
has p = p = 0.001. The initial magnetic field is uniform, 
B = (0.1,0.0,0.0). Figure [g] shows the 2D solution at t = 4 
for rj = 0.018 and rjd — 1/k = 0.18. It exhibits the same 
features as the ideal MHD solution of a similar test problem 
HKomissarov 19 99 ) which is expected given the low value of 
7] and shows nothing that could be suspected as artifacts. For 
more detailed future comparisons with other codes figure [7] 
shows slices of the solution along x = and y = 0. 

The same problem has been used to test the 3D code 
with identical results. 

6.2.2 Strong spherical explosion 

Finally, we tested our 3D code on the problem of spherical 
explosion. All parameters of the explosion are the same as 
in the cylindrical case with exception of the explosion zone 
- now this is a sphere of unit radius. The computational 
domain is (-6.0, +6.0) x (-6.0, +6.0) x (-6.0, +6.0) with 140 
equidistant grid points in each direction. Figures |8] and [9] 
show the numerical solution for 77 = 0.0257 and r^d = l/«; = 
0.257 a.t t = 4.0. The general structure of the solution is 
similar to that of the cylindrical case but with much stronger 
central rarefaction. One qualitatively new feature is the non- 
vanishing electric charge density (top right panel of figH]). 
Given the axial symmetry of the problem one expects the 
solutions to be the same in the planes z = and y — 0. 
Figure [9] shows that this is indeed the case. 



7 CONCLUSIONS 

We have constructed a multidimensional upwind scheme for 
resistive relativistic magnetohydrodynamics. At the moment 
only the case of scalar resistivity has been implemented and 



more work has to be done to incorporate the case of ten- 
sor resistivity. The results of test simulations show that the 
scheme is robust in the regime of small to moderate magne- 
tization, which can be described by the ratio of the electro- 
magnetic energy density to the total mass-energy density of 
matter. The regime of high magnetization is still problem- 
atic as the the truncation errors for the energy-momentum 
of matter become large often making impossible to convert 
the conserved quantities into the primitive ones. This is a 
well know problem of all conservative schemes for relativis- 
tic MHD. Apart from this drawback the scheme can handle 
equally well both resistive current sheets and shock waves 
and thus can be a useful tool for studying phenomena of rel- 
ativistic astrophysics that involve both colliding supersonic 
flows and magnetic reconnection. 
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Figure 9. Strong spherical explosion. In both panels stars show the solution along the Left panel: logiop, gas pressure. Right panel: q, 
electric charge density. 
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